# Establish directories

in.dir = '/home1/99/jc152199/brt/data/'
setwd(in.dir)
out.dir = '/home1/99/jc152199/brt/data/'

# Load the gbm library to perform Boosted Regression Tree Analysis
necessary=c('gbm')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Read in data to model

model.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv',header=T)

# Randomly shuffle model.data before reading into gbm()

model.data = model.data[c(sample(c(1:nrow(model.data)), nrow(model.data),replace=FALSE)),]

# Establish a list of parameters for the BRT model, including learning rate, tree complexity, and training fraction

lr.list = c(.25,.1,.01,.001,.0005) # Learning Rate
tc.list = c(1,3) # Tree Complexity
tf.list = c(.1,.25,.5,.75) # Training Fraction
nt=5000 # Number of Trees

# Establish a dataframe to write summary data into

valid.error.data = data.frame(learning.rate=NA,train.fraction=NA,iteration=NA,optimal.tree.num=NA,min.valid.error=NA)

	for (tc in tc.list)
	
	{

		for (tf in tf.list)
		
		{

			for(lr in lr.list)
		
			{
			
			# Run gbm model
			
			brt.gbm = gbm(formula = micro_max~AWAP_max+coastdist+fpcmean+fpcvar+roaddist+solar, distribution = "gaussian", data = model.data, n.trees = nt, interaction.depth = tc, shrinkage = lr, train.fraction = tf, bag.fraction = 0.5, cv.folds=10)
			
			# Determine optimal tree number using cross-validation on the model data
			
			op.tree.num = gbm.perf(brt.gbm, plot.it = FALSE, oobag.curve = FALSE, overlay = TRUE, method='cv')
			
			# Create a data frame of model summary info

			t.data = data.frame(learning.rate=lr, train.fraction=tf, iteration=which(brt.gbm$valid.error==min(brt.gbm$valid.error)), optimal.tree.num=op.tree.num, min.valid.error=min(brt.gbm$valid.error,na.rm=T))
			
			# Bind to blank data frame
			
			valid.error.data = rbind(valid.error.data,t.data)
			
			}
		
		}
		
	}	

# Write out summary data
	
write.csv(x=valid.error.data, file='valid.error.summary.csv', row.names=F)

#End	
